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TECHNICAL MEMORANDUM 


NUMERICAL STABILITY OF AN EXPLICIT FINITE DIFFERENCE SCHEME FOR 
THE SOLUTION OF TRANSIENT CONDUCTION IN COMPOSITE MEDIA 

INTRODUCTION 


Stability of numerical schemes plays a key role in the computer solution of 
transient heat conduction. A stability analysis relates permissible combinations of 
ffi'id spacing and time step. For a given grid spacing, the time step must be 
smaller tlum a certain value for an explicit numerical scheme to be stable. A time 
step as close to the maximum ns possible is desirable to reduce computer run times 
and cost. 

In composite media, the tivinsient heat conduction equation is not valid at inter- 
faces between media with different heat conduction coefficients. In this case, a heat 
balance technique can be used to drive tm explicit finite difference scheme. The 
derivation is presented herein, and a theoretical stability analysis is performed. 
Because of the complexity of the resulting relationslups , a computer graphics code 
was developed to allow easy determination of an optimum time step. A sample problem 
is examined, and graphics output is presented. 


DERIVATION OF THE FINITE DIFFERENCE SCHEME 

Figure 1 illustrates the general problem that is addressed herein. The following 
assumptions are made: (1) Media interfaces lie along two directions that are orthog- 
onal, e.g. . vortical or horizontal, as shown in Figure 1. (2) Interfaces lie along 
lines of grid points; and when more tlian two media intersect, they do so at a grid 
IK)int. (3) Media interfacial contacts are perfect, i.e. , infinite interfacial conductance, 
(4) Cartesian coordinates are used. (5) Heat conduction is in two dimensions. 

The preceding assumptions were made for simplicity of presentation and are not 
absolute restrictions on the usefulness of the methods described. All results can be 
generalized so that none of the preceding assumptions is required. 

Figure 1 shows a small element centered around grid ix)int (i, j). An explicit 
finite difference ecpiation is desired involving the temperature at (i. j) and surround 
ing grid points at time t and at time t 4 At. where t is the time step. The desired 
result is accomplished using a heat balance. The heat balance is in the form 


Chanifv of iriterria/ energy within the element during the time step - 
heat eonduction into the element during the time step. 


A heat balance form is used because the heat conduction equation is not valid at the 
interface between media. In equation form, the preceding relation is 



MEDIA 

INTERFACE 



Figure 1. Finite difference grid arrangement for the general case 

of four intersecting media. 
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where 

T - temperature 

i = subscript corresponding to x direction 

j = subscript corresponding to y direction 

n = superscript corresponding to the time step 



Ax - X grid spacing 
Ay = y grid spacing 
At = time step 

kj, kg, kg, = thermal conductivity of the media in quadrants 1 through 4 
Pl, pg, Pg, P 4 = density of the media. 

The terms in equation (1) are 


Change in the internal energy of the element during time At = heat 
flow into the element from the left-hand side during time At + heat 
flow into the element from the bottom + heat flow into the right-hand 
side + reat flow in from the top. 


Equation (1) is an explicit form; i.e. , the temperature one time step ahead can 
be determined from temperatures at the grid point and adjacent points at the current 
time step. To see this more clearly, equation (1) is rewritten in the following form: 
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If the media in the four quadrants are the same, equation (2) reduces to 
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Observe that equation (3) is just the finite difference form of the heat conduction 
equation. The assumption of steady state heat flow with Ax = Ay reduces equation 
(3) to 



^ Vl.j '^i>i-l 


(4) 


The superscript is dropped from equation (4) because it is superfluous. Observe 
that the preceding equation is a well-known finite difference form of LaPlace's equa- 
tion, i.e. , the steady state heat flow equation. 

The finite difference form for the heat balance equation was derived and is 
given by equation (2), This equation is in explicit form and, hence, is only condi- 
tionally stable. An inequality relating At, Ax, Ay, and the material parameters must 
be found so that the marching technique suggested by equation (2) is stable. For 
equation (3), the relationship can be shown to be [1] 


k At ^ 1 

2^4 

cp(Ax)^ 


(5) 


The normal procedure for obtaining stable solutions is to select Ax for adequate reso- 
lution of the temperature field and use equation (5) to calculate the maximum At. The 
relationship analogous to equation (5) for the heat balance equation will be derived in 
the next section. 


STABILITY OF THE EXPLICIT FINITE DIFFERENCE EQUATION 


The stability analysis of the finite difference equation (2) will follow the pro- 
cedure outlined by Richtmyer and Morton [1]. First consider the real solution of 

equation (2) given by T.”. . The actual solution calculated on the computer will be 

n 

where i''. . is a small error that may be due to roundoff or trunca- 
* •] 

tion, function evaluation error, etc. 

tion (2), it can be seen that (S is a solution to the same equation as T.” since 

^>1 i»j 

Reference 1 shows that the error 


f = T." + (S." 

1 .] 1.1 »0 


n 


If the solution with error is plugged into equa- 


T. . is a solution of (2) and hence will cancel out. 
* d n 

function • can be represented in the form 
* 0 


r n ^ n i(i3Ax + JyAy) 
i.j ^ 


( 6 ) 


where i = / - 1 
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Here C is, in general, a complex amplitude factor; and means ^ raised to the 
nth power (not ^ superscript n). 3 and y can be considered free parameters that 

can change independently of each other. Substituting equation (6) into equation (2) 
yields, a /ter some rearrangement. 


r ^ ^ ^4^ "I 

^ = 1 + a I - — (cos B + cos G - 2n 

+ 1 a 


^4 ‘ ^ ^2 ■ *^3 ■ ^4 

S sin B + =: 


) “I 

- sin G , 


(7) 


where 


^ ^ 4At 

AxAy(Cjpj + C2P2 + CgPj + c^p^) 

B = 3Ax 

G = yAy . 


In the example of the discretized form of the heat conduction equation, C is real. 

For the current form, C is. in general, complex. For a given grid point, the ampli- 
fication factor for any error is given by equation (7). If the modulus of ^ is greater 
than unity, any error present can amplify and disrupt the numerical solution of ♦’ 3 
temperature field. To assure stability of the solution, 1 must be satisfied at 

every grid point in the interior of the composite region and at every boundary point 
(if heat flux or radiation boundary conditions are specified). 

For a given g^rid point in a composite medium, C is a two-parameter function of 
B and G in the complex plane. If Ax and At are selected and 1 1 for every 

value of B and G, then the stability condition is not violated. 

To better understand equation (7), fix G at some arbitrary value. Then the 
following representations can be made: 


Re[ r] = Cj cos B + C 2 

Im[^] = Cg sin B + c^ , (8) 

where 

RelfJ = real part of f, 

Im[fJ = imaginary part of r. 
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Then 
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c 
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= sin B 


Finally , 


(X - C2>‘ 


(y - Cd>‘ 


( 10 ) 


This is the equation of an ellipse centered in the complex ^ plane at ^ = C 2 + i 
with semimajor axes c^ and Cg. Reference to the definitions of and Cg indicates 

that the semimajor axes are constants for a given composite medium; i.e., they are 

not functions of G. c« and c- are functions of G and can be written in the form 

2 4 


ir 


Cg = COS G + 

c^ = sin G 


( 11 ) 
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where 


a(kj * kj - “3 ■ ’'4* 
= 2 


Equation (11) is the parametric form for an equation of an ellipse. Equations (10) 
and (11) represent a series of ellipses centered on points on an eUipse. At this 
point, two approaches can be taken. The first to be considered involves trying to 
find an envelope for the series of ellipses defined by the preceding two equations. 
Determining an equation for the envelope involves combining equations (10) and (11) 
in the following form: 


2 2 
(x - a, cos G - a_) (y - 3, sin G) 

f(x.y.G) = 2 — 2 ' ^ ‘ 

*^1 °3 


If an envelope exists, it is necessary that it satisfy the following two equations (see 
Reference 2, for example): 


f(x,y.G) = 0 

f^(x.y.G) - 0 . (13) 


The subscript G refers to partial differentiation. Conditions (13) become 


2 2 
(X Uj cos G - 1X2) (y - 3j sin G) 



tan G + kj sin G - k 2 = 0 


(14) 


where 
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The usual procedure for solving for the envelope involves elimination of G from the 
two equations in (14). Because of the transcendental nature of these equations, 
elimination of G is difficult. Rather than proceed in an attempt to obtain an analytical 
solution, a computer program was written to plot up the series of ellipses. This pro- 
gram is described in the following section. 


COMPUTER GRAPHICS STABILITY DETERMINATION 


The graphics program is based on equations (8) and (11), which are presented 
here in slightly different form for convenience, 

X = Cj^ cos B + cos G + 

(15) 

y = Cjj sin B + sin G 


This represents a two-parameter family of curves. The finite difference form of 
equation (2) is stable, if the following inequality is true: 


= / 




‘-2 . 

y i 


(16) 


The preceding two -parameter family is plotted by holding G constant and varying B 
between 0 and 2n. This plots one ellipse. After G is incremented between 0 and 2n, 
a sedes of ellipses is drawn that allows visual determination of the C envelope. The 
program listing in the appendix draws the circle Ic | = 1 and the ellipses. If a part 
of the envelope falls outside the unit circle, instabilities can be expected. If this 
occurs. At must be reduced and a new graph drawn. When the envelope just stays 
in the unit circle, the corresponding value of At is nearly optimum. 

To better illustrate the method, the following section outlines an example prob- 
lem and indicates how optimum time increments can be obtained. 
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EXAMPLE STABILITY ANALYSIS 


The example ^metry of Figure 2 is used to illustrate principles outlined in the 
preceding sections. The composite medium is made up of copper, plexiglass, cork for 


insulation, and Dow -Corning 200 Series oil (viscosity 1000 centistokes). Although the ! 

oil is a fluid and subject to convection, its high viscosity and the low temperature f 

gradients of the problem (low Rayleigh number) mean that most heat transfer is b- I 

conduction. Ax and Ay are 0.1 cm for the example. Given these values for t.y c uu 
Ay, the object of the stability analysis is to determine the maximum At for a st 'le i 

solution . 




Figure 2. Composite medium geometry. 

Each grid point can be classified into one of several categories, as indicated 
in Figure 3. To assure stability, each of these categories must be tested. For the 
example problem, stability considerations ere dominated by copper because of its high 
conductivity relative to the other three materials of the composite medium. Exper- 
ience shows that the only categories of interest are those containing copper in at 
least one of the four quadrants. 

A trial and error procedure indicates that At - 0.002 sec is nearly the maximum 
or optimum value. Figure 4 sh^ws the stability plots for the interiors of each of the 
four materials. C is real :n . e interiors and becomes significant only for the case of 
copper. The plots of Figure 3 show the remaining categories listed in Figure 3. 
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Figure 3. Classes of grid point arrangements that must be 
considered for stability. 

Because of the dominance of copper, all curves for the interfaces appear to 
circles with centers on the real axis. Although not belonging to one of the cate- 
gories, the plot of Figure 6 is included to slx)w more fully the general character of 
the stabMity equations. Even this plot appears to be composed of circles. Again, 
the dominance of copper is evident. If two dominant materials existed with compar- 
able but unequal conductivity, a series of ellipses could be expected. 

To check the results of this stability analysis, the thermal code based on 
equation (2) was programmed. For a value of At = 0.002 sec, the code was stable. 
For At = 0.00225 sec, instabilities soon developed in the solution. 

If inequality (5) for copper is used to find the maximum value of At, the 
result is At = 0.0021665 sec. Figure 7 shows the copper category stability plots. 

The amplification factor closely approaches the unit circle. Figure 7a serves as a 
check on the code. These findings suggest that stability is determined by inequality 
(5) applied to the dominant material of the composite medium. The restilt is by no 
means proved, however. 
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RE (XI) 
5(«) 



2 - 1 . 0 . 


RE (XI) 
S(f) 


I'ipure 5 


TIME INC ° 
.00200 SEC 



1 . 2 . 


MATERIALS 
1*CU 
2=C0RK 
3==DC OIL 
4=PLEX 


TIME INC • 
.00200 SEC 


(Concluded) . 







garc 6, Hypothetical p*id point with a different material 
in cacli of the four quadrants. 




Fijjure 7, Stability plots for the optimum At. 














MATERIALS 

1- CU 

2- CORK 

3- OC OIL 

4- PLEX 


TIME INC > 
M217SEC 


2. 1. P- 1. 2- 

RE (XI) 

7|i) 

Figure 7. (Concluded). 




SUMMARY 


A computer graphics technique was derived that was useful for determining the 
maximum time step in an explicit finite difference equation describing transient con- 
duction in a composite medium. An example was used to illustrate the method. For 
the example, two codes were programmed. The first was the explicit, conditionally 
stable thermal code based on equation (2). The second, or stability code, was a 
computer graphics program for determining the optimum time increment which was an 
input for the thermal code. Values of At for which the stability code indicated 
stability (At = 0.002 sec) and instability (At = 0.00225 sec) were input to the thermal 
code. The results from the thermal code were consistent with the predictions from 
the stability code. 

For the example problem, an accurate stability prediction could be achieved by 
applying inequality (5) to the dominant materitil, which in this case was copper. 

Tlhs result may or may not be a universal one. 

The stability of the conditionally stable explicit finite difference equation (2) 
was controlled by the dominant material in the example problem. For the other three 
materials, a much larger maximum time step would be calculated from (5). Two 
possibilities exist for overcoming this problem. One is to use a larger Ax and Ay in 
the dominant material. For the example, this procedure was not feasible because 
sufficient resolution could not be achieved. The second possibility is to use a 
different time step in the different materials. This possibility means that, while the 
solution is updated every 0.002 sec in the copper of the example, it is only updated 
every 0.2 sec in the other materials. This variable time method should be explored. 
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APPENDIX 


This program listing is written in Hewlett-Packard (HP) Fortran 4X, which is 
described in Reference 3. Except for the multiple statement lines in which statements 
are separated by the dollar sign ($), the program uses standard Fortran IV. The 
graphics subroutines are part of a special package called Graphics 1000, which is 
described in Reference 4. 
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0030 
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0034 
0 035 
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0038 
0038 
004 0 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 
0 048 

0050 

0051 
0 052 
0 053 
0054 
0 055 
0 056 
0 057 
0 058 


FTN4 ,L 

PROGRAM TSThP 

C iti-*i*4‘*i4i*iti‘ti**'ti**** 4<>f ff 1-* ******* **■*■** * ** 4i*****-**H'*****i4iH'-ti't'****'¥i¥*‘t"*‘t>*>¥ 


C ♦><<!♦ ♦♦♦ 

c -f** PROGRAMMER: WARREN C AMPBELL/3- 1 886 ■*♦* 

C 

C *++ PROGRAM PURPOSE; TO CALCULATE STABILITY FOR FINITE DIFFERENCE *** 

C ■»•*•+ SCHEME FOR SOLUTION OF TRANSIENT CONt-UCTION EClUATION IN 

C COMPOSITE MEDIA WHOSE INTERFACES ARE VERTICAL OR HORIZONTAL. 

C i< * ★ 

I*' .«i « il 4l K 

C *f*- DEFINITION OF VARIABLES- *** 

C **» 

C NOTE; ALL UNIT? hRE CGS **>• 

C ^ ♦ 

C *++ RriOCU =• DENSITY OF COPPER ■*** 

C *■** RHOCK “ DENSITY OF CORK ♦♦♦ 

C RHODC = DENSITY OF DOW CORNING 200 SERIES OIL *** 

C RHOPL » DENSITY OF PLEXIGLASS *** 

c + + ♦ ecu = SPECIFIC HEhT OF COPPEP >»•♦>*■ 

C CCK = SPECIFIC HEAT OF CORK *** 

C «♦"*• CDC = SPECIFIC HEAT OF DOW CORNING 200 SERIES OIL 

C CPL « SPECIFIC HEAT OF PLEXIGLASS 

C KCLI = THERMAL CONDUCTIVITY OF COPPER 

C *■*"*■ KCK = THERMAl CONDUf T I V I TV OF CORK ♦♦♦ 

C kDC »= THERMAL CONDUCTIVITY OF DOW CORNING 200 SERIES OIL 

C *■*■->■ KPL = THERMAL CONDUCTIVITY OF PLEXIGLASS 

C ♦'♦■I RH01 ,RH02.RH0>.RH04 = DENSITY IN 1 ST , 2ND . 3RD , 4TH QUADRANT t** 

C C1,C2,C3,C4 = SPECIFIC HEAT IN " " " " " h..*.* 

L k1,K2,K3.K4 = thermal CONDUCTIVITr If) iST ETC LJUmDRANT *** 

C *+♦ DELX = X INCREMENT IN FINITE DIFFEFfENCE EQUATIONS >*■♦>• 

C DELY = Y '* " " " " 

C *** DELT = TIMF INlRFHFNT >♦♦=* 

C A - DEFINED PARAMETER 

C XI = COMPLEX amplification FACTOR **■* 

c X = REAL Part of the amplification Factor -k** 

c ♦ Y = imaginary part ok amplification factor 


DIMENSION IPLTR.; 182 ).Cv4 ),RHO':4>. IM' 4 > 

REAL K'. 4 ) , K C U , K C k . KDC . KPL 

P 1=3. 14158 

DELX=0 . 1 

DELY=0 . I 

WRITE', I , 888 ' 

888 FORMAT. " ENTER TINE INCREMENT i! " < 
PEaIk 1 ) DELT 

C ♦ ♦ COPPER PARAntTER:- 
KCU=0 ,8512 
RHOC U = S . 86 
CCU=0 , ''^2 

C •* f- * CORK PARaMFTcRS *t* 

K C K = u 0 i< 0 1 0 i 
RHOc 0,1 6 02 
CCK=0 04 

C >’*> Dl 200 OIL PnPh.:iE rEP'I- *** 
fDC i'lOOT” 
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0059 

0060 
0061 
0062 

0063 

0064 

0065 

0066 
006? 
0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 
0073 

0079 

0080 
0081 
0082 

0083 

0084 

0085 

0086 
0087 
0038 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0 093 
0099 

01 00 
01 01 
01 02 
01 03 
01 04 
01 05 
01 06 
01 07 
01 08 
01 09 
01 1 0 
01 1 1 
0 1 1 i 

0113 

0114 

0115 

0116 
0117 
01 18 


RHODC-0.971 

COC-0.379 

C ENTER PLEXIGLASS PhRhMETERS ♦♦♦ 

KPL-4 .475E-4 
RHOPL»1 . 19 
CPL-0.35 

C *** ENTER MATERIAL FOR EACH OF THE FOUR QUADRANTS 
UR1T£< 1 , 1000> 

1 000 FORMAT< " CU«1, C0RK«<2. DC 0IL«3, PLEX«4, QUAD 1 MAT-<1 "> 

READ< 1 ,<•> IM< 1 > 

URITE< 1 , 1 1 00> 

1 1 00 FORMAT< " QUAD 2 MAT»^ '• > 

READ<1,*) IM<2) 

URITE< 1 , 1200> 

1200 FORMATC QUAD 3 MAT»4i " 

READ< 1 ,♦) IM<3) 

WRITE< 1 , 1300> 

1300 FORMAT< " QUAD 4 MAI ' > 

REAO<1,>*) IM<4> 

DO 130 I»1,4 

IF< IM< I >.NE, 1 ) GO TO 1 00 
K< I >»KCU 
RHO< I )»RHOCU 
C< I >=CCU 
GO TO 130 

100 IF< IM< I ).NE.2 ' GO TO 110 
K< I )=KCK 
RHO< I >=*RHOC.K 
C< I >*CCK 
GO TO 130 

110 IF< IM< I >.NE.3:) GO TO 120 
K< I 

RHO< I ' = RHODC 
C< I :»=CDC 
GO TO 130 
120 K< I )=I-'PL 

RHO< I )*RHOPL 
C< I >=CPL 
130 CONTINUE 

m» 4 . ♦DELT/'i DELX *DEL Yt". C’. 1 »*RHO< 1 > ♦•C', 2 > t-RHOt 2 J + ?• >+RHOi; 3 }*C\ 4 >■•• 

1 RHOs 4 > > .) 

C f** SET UP PLOTTER 
WRITE-: 1 , 1 1 1 1 > 

1111 FORMAT-'" ENTFR ID NUMBER > 

READ< 1 , * ) ID 

WRITE< 1 , 1 1 12 > 

1112 FORMAT-:" ENTER LU NUMBER ' 
read* 1 , ♦ ) LU 

CALL PLOTR-: IPLTR . ID, I , LU •> 

CALL VIEWP< IPLTR, 15. ,85. . 15. ,85. - 
CALL UHNDW-; IPLTR, -2 , 2 . , -2 . , 2 . 

IPEN=1 

CALL PEN< IPLTR, IPEN > 

C ALL L AXES'. I PL TP , - , 2 . 0 . 2 . 0 . , 0 . , 5 . . 5 . , I . > 

CALL FRAME< IPLTP > 

CALL MOVE-. IPLTR, 1 0. > 

DO 200 1=1,101 
hRG*'. 1-1 :^2 . *P I / 1 00 . 

X*lOS-. ARG ) 
y«SlN'. ARG ) 
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0119 


CALL DRAU< IPLTR,X, «■ - 

0120 

200 

CONTINUE 

0121 


CALL PENUP(IPLTR) 

0122 


DO 210 I ■ 1 , 2 1 

0123 


CALL PEN< IPLTR, IPEN ) 

0124 

' 

G-<1-1 . )*2.*PI/20. 

0129 


DO 220 J-1.51 

0126 


B«< J-1 )-*2.*PI/50. 

0127 


X"1 . +A*< 0 . 5*( K< 1 >+«< 2 )+K<! 3 >+K< 4 > )>*^ C0S( 8 >+COS< G >-2 . ) > 

0128 


Y«Ai*<<K< 1 )+K<4)-K<2>-Kc3^)'*.5i-SlN>:B)+<K-: 1 2 i-K< 3 4 ) )* . 5 

0129 


1-*8IN<G>> 

0130 


CALL DRAU< IPLTR,X,V ) 

0131 

220 

CONTINUE 

0132 


CALL PENUPdPLTR) 

0133 


IPEN-IPEN-H 

0134 


IF< IPEN.LE.4) GO TO 210 

0135 


IPEN«1 

0136 

21 0 

CONTINUE 

0137 


CALL VIEWP'; IPLTR, 0. , 152. . 0. , 1 00, ; 

0136 


CALL WINDW". IPLTft, 0, , 152. , 0. , 1 00. ; 

0139 


IPEN«1 

0140 


CALL PEN< IPLTR, IPEN > 

0141 


CALL MOVE< IPLTR,50. , 10. ) 

0142 


CALL CPLOT-; IPLTR.-3, 0. > 

0143 


CALL LABEUIPLTR) 

0144 


WRITE<LU,4000> 

0145 

4000 

FORHAT< "9E<XI )" ) 

0146 


CALL MOVE< IPLTR, 1 0. ,50. .» 

0147 


CALL CPLOT'. IPLTR, 0. ,-3. > 

0148 


CALL LOIR-; IPLTR, 1,57) 

0149 


CALL LABEL< IPLTR) 

0150 


WRITE<LU,41 00 ) 

0151 

4100 

FORMAT< "IM-;XI )“ .) 

0152 


CALL LDIR-; IPLTR, 0. > 

0153 


CALL MOVE-; IPLTR, 50 90 . ) 

0154 


CALL CPLOT' IPLTR, -13. , 0, .) 

0155 


C«LL LABELC IPLTR; 

0156 


WRITE<LU,42 0 0.) 

0157 

42 00 

FORMATC "AMPLIFIChTIuN FhlTOR PluT" ■ 

0158 


CALL MOVE< IPLTR, 90 . , 85 . ' 

0 1 59 


CALL LABELC IPLTR; 

0160 


WRITE< LU, 43 0 0 ) 

0161 

4300 

FORMAT^ “MHTERI«LS“ ) 

0162 


WRITE< LU, 44 0 0 ) 

0163 

4400 

FORMATS" 1*CU" ) 

0164 


WR1TE< LU, 4500 '> 

0165 

45 00 

FORMAT-;" 2= CORK") 

0166 


WRITE<LU,4600 ) 

0167 

46 11 0 

FORMATS" 3-=C’C OIL") 

0168 


WRITE< LU, 4700 > 

0169 

4 7 0 0 

FOPMATc^ " 4 = PLEX") 

017 


CALL MOVE'. IPLTR, 1 1 0. ,75. / 

0171 


CALL DRAW-' IPLTR, 1. ;0., 75 .) 

0 1 72 


Call MOVE*. IPlTR , 1 1 5 . , 70 , 

0173 


CALL DRAW-: IPLTR, 1 15, . 80 . ! 

0174 


CALL MOVE', IPLTR, 1 17, .77, ) 

0 1 75 


CALL LABEL-. IPLTR ; 

0176 


WRITE< LU, 480 0 > IM-: 1 > 

0177 

48 00 

FORMATC I 1 ) 

0178 


CALL MOVE', IPlTR, 112 .77, ) 



0179 

CALL LABEL’; IPLTR/ 


0180 

WRITE<LU,4800> 


Oidi 

CALL MOVE< IPkTR, 1 l2 

, 72 , 

0182 

CALL LABEL<IPLTR> 


0183 

UIRITE<LU,480 0> IM':3 > 


0184 

CALL MOVE< IPLTR, 117. 

,72. ) 

0185 

CALL LABEL< IPLTR) 


0186 

I<IRITE<LU< 480 0) IM’;4) 


0187 

CALL MOVE< IPLTR, 90. , 

30. ) 

0188 

CALL LABEL< IPLTR) 


0189 

WRITE<LU,4900 ) DELT 


0190 

4900 FORMATCTIME INC = " 

F8.5“ 

0191 

CALL PEN< IPLTR, 0) 


0192 

STOP 


0193 

END 


0194 

END* 



(f 
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